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ABSTRACT 

A  theory  combining  nonlinear  continuum  elasticity,  inelasticity,  thennodynamics, 
damage  mechanics,  and  fragmentation  is  fonnulated.  The  model  is  applied  to  study  concrete 
subjected  to  high  rate  loading  as  occurs  during  ballistic  impact.  Two  thermodynamically 
motivated  methods  are  postulated  within  this  theoretical  framework  for  quantitatively 
characterizing  the  mass  and  velocity  distributions  of  the  post-impact  debris  field,  one  based 
upon  a  local  energy  balance  and  a  second  following  global  entropy  maximization.  Here  the 
concrete,  a  composite  mixture  of  mortar  and  granite  aggregate,  is  regarded  as  a  homogeneous 
continuum  prior  to  fragmentation.  However,  the  composite  nature  of  the  microstructure 
directly  influences  model  parameters  dictating  the  mean  fragment  dimension,  here 
specifically  related  to  the  coarse  aggregate  size.  Standard  continuum  elements  represent  the 
intact  solid  and  particles  describe  eroded  material  in  numerical  implementation  of  the  model. 
The  impact  of  a  metal  sphere  on  a  thin  concrete  target,  and  the  subsequent  motion  of  the 
resulting  cloud  of  concrete  fragment  debris,  are  simulated.  Fragment  size,  speed,  and  kinetic 
energy  statistics  predicted  by  the  two  methods  are  compared. 

Keywords:  damage  mechanics,  thermodynamics,  fragmentation,  concrete. 

1.  INTRODUCTION 

Urban  structural  materials  such  as  concrete,  mortar,  and  cinder  block  may  undergo 
complex  set  of  deformation  mechanisms  when  subjected  to  impact  loading.  Since  they 
nominally  contain  significant  initial  porosity,  for  example  on  the  order  of  10-20%  in 
standardized  concretes  [1,2],  these  materials  are  labeled  here  as  'crushable'.  The  goal  of  the 
present  study  is  development  and  implementation  of  a  kinematically  and  thermodynamically 
self-consistent  theory  accounting  for  coupled  deformation,  damage,  and  fragmentation 
mechanisms,  specifically  amenable  to  brittle  crushable  solids.  In  such  materials,  initial 
porosity  induces  pressure  dependence  in  the  effective  bulk  modulus,  with  the  stiffness  of  the 
material  increasing  as  the  pores  are  compacted.  Additional  difficulty  in  precisely  modeling 
concrete  behavior  follows  from  the  variability  in  mechanical  properties  such  as  flow  stress 
and  fracture  toughness  with  microstructure  constituents  [3],  processing  conditions,  and  age. 

Previous  modeling  approaches  most  relevant  to  the  present  work  are  mentioned  here. 
Hohnquist  et  al.  [4]  developed  a  constitutive  model,  with  the  pressure-volume  response 
following  data  from  [1],  and  a  plasticity  and  damage  framework  similar,  yet  not  identical,  to 
one  used  previously  for  metallic  materials  [5],  with  failure  criteria  based  on  cumulative  strain 
measures  [6].  The  numerical  approach  followed  here  is  that  of  Johnson  et  al.  [7,  8],  whereby 
continuum  finite  elements  are  converted  via  a  Generalized  Particle  Algorithm  (GPA)  to 
particle  nodes  when  specified  erosion  criteria  are  met.  However,  infonnation  regarding  sizes 
of  individual  fragments  has  not  historically  been  available  from  Smooth  Particle 
Hydrodynamics  (SPH)  or  GPA,  since  the  mass  of  each  particle  is  simply  the  nodal  mass, 
which  in  turn  depends  upon  the  mesh  density.  The  issue  is  resolved  here  by  incorporating 
thennodynamics-based  fragmentation  theory  into  the  constitutive  framework  enabling 
calculation  of  non-uniform  mass  distributions. 

Two  methods  are  derived  in  the  present  work  for  calculating  fragment  size  and 
velocity  distributions,  both  compatible  with  the  laws  of  thermodynamics  and  momentum 
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conservation.  The  first  follows  from  on  a  local  energy  balance  similar  to  that  of  Grady  [9] 
and  Johnson  &  Cook  [10],  but  newly  applied  to  the  class  of  concrete-like  materials  in  a 
manner  more  consistent  with  the  bulk  constitutive  model.  The  distribution  of  fragment 
masses  is  found  in  practice  via  consideration  of  an  energy  balance  with  regards  to  the 
deformation  history  of  converted  particles  in  the  simulation,  with  the  velocity  distribution  of 
the  fragments  associated  with  that  of  the  parent  particles.  In  the  second  approach,  a  theory 
based  on  entropy  maximization  and  classical  statistical  physics  is  implemented  [11,  12],  In 
this  approach,  a  new  joint  probability  distribution  function  for  fragment  mass  and  velocity  is 
derived  consistent  with  energy  and  linear  momentum  conservation 

The  following  notation  applies.  Cartesian  coordinates  are  used  throughout,  with 
summation  implied  over  repeated  indices.  Vectors  and  tensor  quantities  are  represented  with 
boldface  type,  while  scalars  and  individual  components  of  tensors  are  written  in  italics. 

Juxtaposition  implies  summation  over  two  repeated  adjacent  indices  (e.g.  (AB)^  =  AacBcb). 

The  scalar  product  of  vectors  is  represented  by  (e.g.  a*b  =  a“ba).  The  colon  denotes 

contraction  over  repeated  pairs  of  indices  (e.g.  A  :  B  =  tr(ArB)  =  AabBab ,  where  V  is  the 

trace  operation,  and  C :  A  =  Cahcd Acd  ).  Superposed  -1,  T,  and  denote  inverse,  transpose, 
and  material  time  derivative,  respectively. 

2.  CONTINUUM  THERMOMECHANICS 

The  theoretical  framework  postulated  here  differs  from  many  existing  models  in  its 
implementation  of  a  multiplicative  decomposition  of  the  deformation  gradient  into  elastic  and 
inelastic  components  [13]  and  its  strict  consistency  with  the  balance  laws  of  continuum 
mechanics  and  thermodynamics.  More  specifically,  thennodynamic  relations  and  evolution 
equations  for  porosity  and  damage  are  formulated  in  the  spirit  of  Coleman  &  Noll  [14]  and 
Coleman  &  Gurtin  [15]. 

The  defonnation  gradient  F  is  split  multiplicatively  as 

F=  dx/dX  =  FeFd  ,  FaA  =  9xa / dXA  =  FEaaFDaA ,  (1) 

where  x  and  X  are  spatial  and  reference  coordinates  in  three  dimensional  Cartesian  space,  FE 
is  the  recoverable  elastic  deformation,  and  F°  is  the  irreversible  deformation  associated  with 
defects  such  as  micro-cracks  and  voids  within  the  material.  Deformation  maps  and 
corresponding  material  configurations  corresponding  to  Eq.  (1)  are  depicted  in  Fig.  1.  The 
velocity  gradient  L  follows  as 


L  =  dx/dx  =  FeFe-1/  +  FeFdFd-1Fe  1  .  (2) 

le  ld 

Irreversible  volumetric  deformation  associated  with  pore  collapse  in  crushable  materials  is 
described  by 

(p  =  JD~l- 1,  JD  =  detFD ,  (3) 

where  (p  is  the  volume  reduction  upon  crushing,  positive  when  the  volume  is  reduced.  The 
total  inelastic  velocity  gradient  from  (2)  then  becomes 

Ld  =  FEFDFD1FE1  =  Ld  -  q>{\  +  <p)~l  1 ,  (4) 
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where  LD  is  the  deviatoric  inelastic  velocity  gradient  in  spatial  configuration  B  and  1  is  the 
identity  tensor.  An  elastic  strain  tensor  and  a  scalar  measure  of  volumetric  elastic  strain  in 
the  intermediate  configuration  B  of  Fig.  1  are  defined,  respectively,  as 

2EE=FErFE-l,  &E  =  trEE .  (5) 

Standard  balances  of  mass,  linear  and  angular  momentum  are  written  in  localized  fonn  as 

p0  =  pJ  —  p  det  F ,  diva  +  f  =  px,  o  =  or,  (6) 

where  f  is  the  body  force  vector  per  unit  spatial  volume,  o  is  the  Cauchy  stress  tensor,  and 
div  denotes  divergence  in  the  spatial  frame.  The  local  spatial  energy  balance  is 


pe  =  o  :  F  -  div q  ,  (7) 

with  e  the  internal  energy  per  unit  mass  and  q  the  heat  flux.  The  local  second  law  of 
thennodynamics  is  written 

o:\^-dl(\»dJ)>p[y  +  dT^,  (8) 

where  6  is  the  absolute  temperature,  r]  is  the  specific  entropy  per  unit  mass,  and  the 
Helmholtz  free  energy  is  y/  =  e-drj .  The  spatial  gradient  operator  is  denoted  by  dx . 

The  Helmholtz  free  energy,  on  a  per  unit  mass  basis,  exhibits  the  following  functional 

form: 

yy  =  yy{^,(p,e,D),  (9) 

where  D  is  a  scalar  internal  state  variable  representing  cumulative  damage.  Stress-strain  and 
temperature-entropy  relations  are  then  deduced  as 


o  =  F> 


e  „  dy/  ET 


5E1 


v  =  -p 


dy / 
~d0 


(10) 


The  dissipation  inequality  (8)  is  then  given  by 


o  :F°  + 


p  dy/ 
(l+</9  P  d(p  ) 


dp  —  p—^—D  >  0  'q *8  6 
QD  x 


(11) 


where  the  Cauchy  pressure  3 p  =  -tro .  Assuming  isotropic  heat  conduction  in  the  spatial 
frame,  i.e.  q  =  -kBj9  ,  with  k  the  scalar  thennal  conductivity,  and  defining  the  specific  heat 
parameter  c  =  de i 80  =  -6d2y/  /  dO2 ,  the  energy  balance  can  be  written  as 


pcO  -  o  :  Fd  -  p 


dy  e  gV 

dtp  dOdcp 


cp- 


dy/ 

dD 


■0 


d2y/ 

86dD 


D  +  0- 


dy/ 

dm} 


■:  E 


+  div{iddj)') .  (12) 


In  the  present  class  of  urban  structural  materials  (e.g.  concrete,  mortar)  Fd  is 
associated  with  micro-crack  opening  and  sliding,  as  well  as  pore  collapse  during 
compression.  Let  co  represent  the  cumulative  local  micro-cracked  area  per  unit  intennediate 
volume.  Further,  let  D  =  toltoc,  where  <x>c  is  a  material  parameter  denoting  the  maximum 


3 


sustainable  crack  density,  and  constrained  by  the  restriction  0  <  D  <  1 .  The  specific  free 
energy  density  is  more  specifically  formulated  as 

py/  =  K(3E,(p)3l  +  G(\-D)t*  +r(D)  +  Y(d),  (13) 

where  p  =  pJE ,  K  is  the  effective  bulk  modulus,  G  is  the  initial  shear  modulus, 
Ee  =  Ee  -(i9£  /  3 )  1  is  the  elastic  strain  deviator,  r  accounts  for  energy  of  micro-cracks,  and 
Y  describes  the  specific  heat  content.  In  particular,  the  effective  bulk  modulus  takes  the  form 


K  =  -Kc 


^  +  [\kAki9e+\k,A^, 

<Pl  )  \ 2  3  4  )<Pl 


where  KE  is  the  initial  elastic  bulk  modulus,  cpL  is  the  maximum  porosity  reduction  due  to 
compressive  pressure,  and  Kx,  K2,  and  K,  determine  the  pressure-volume  relationship  for 
the  fully  dense  material  at  locking,  when  <p  =  <pL.  The  microcrack  energy  per  unit  volume  is 


-r  =  yco  =  (K2c(ocl2KE)D,  (15) 

where  y  =  K2  /  2K h:  is  the  surface  energy  of  fracture  [16],  with  Kc  the  effective  fracture 

toughness.  Note  that  constant  material  properties  are  used  in  (15)  as  a  first  approximation. 
The  intermediate  second  Piola  stress  is  given,  following  from  (10),  by 

S  =  p-^-  =  JEFE1oF* -t.  (16) 

y  5Ee 

Then  from  (13)  and  (14),  respectively,  hydrostatic  and  deviatoric  parts  of  S  are 

P  =  -trS  /  3  =  -Ke3e  (<pL  -  <p)  /  cpL  -  (K,9E  +  K292E  +  K3$ ) cp  /  <pL ,  (17) 

S  =  S  +  jpl  =  2(l-D)GEE.  (18) 

The  Cauchy  pressure  is  -3 p  =  trc  =  JE^FEaaSaP FE a/3 . 

Deviatoric  plasticity  is  controlled  by  the  flow  potential  Q ,  equated  here  with  the 
scalar  effective  deviatoric  stress  a  =  J(3/Tj6  :  a  : 

LD=i||,  0  =  \A(\-D)  +  B(pla„fJ\  +  C\n{Als„)]<ja,  (p/<70>-f),  (19) 

where  3 A2  =  2L11 :  L°  for  stresses  exceeding  the  elastic  limit,  a  is  the  deviatoric  stress  tensor. 
Material  parameters  from  ref.  [4]  include  A,  B ,  C ,  N ,  s0,  and  a0 ,  and  T  is  defined  as 
-p/<jQ,  with  p  the  tensile  pressure  at  failure.  An  isotropic  inelastic  response  is  assumed 

such  that  the  inelastic  spin  may  be  neglected,  meaning  LD  =  LDr .  Note  also  that 
Ld  «  FDFD1  when  elastic  shape  changes  are  small.  A  2D  slice  of  the  yield  surface  from  (20) 
is  shown  in  Fig.  2  for  the  particular  case  where  0  =  a0 ,  where  al  and  a2  are  principle 

Cauchy  stresses,  and  the  third  principle  stress  cr3  =  0  . 
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The  following  kinetic  equations  dictate  porosity  and  damage  evolution: 


jO  (p<pc;p>pL), 
\d(p)/a0  ( Pc<P<Pl )> 


D  =  KA(l-7TD(p/<Joy}  +  T~l  (~p/<T0  ) ,  (20) 


where  a  ,  k  ,  and  n D  are  positive  constants,  and  the  notation  2  (x)  =  x  +  |x| .  The  pressure  at 
which  inelastic  crushing  commences  is  labeled  pc ,  and  the  locking  pressure  corresponding  to 
<pL  is  denoted  as  p,  .  Note  that  both  <p  and  D  are  always  positive,  i.e.,  irreversible. 

The  model  is  next  demonstrated  to  be  thermodynamically  consistent.  From  ( 1 3)-(  1 5), 


dif / 


E-  1/-tE  .  fiE  .  tE-  1  KcO)c 


-p^  =  JtlGE^:^+J 

dD  2  K 


~pJp=  jE~'  [PoV  +  KA  l2(pL-{KJ2  +  K23Ell  +  K,&2El4)&2El(pL 


Returning  to  (1 1),  enforce  the  stronger  requirements 


„:L °-p^D  + 

dD 


'  p  dy/^ 

— - p—1— 

U  +  P  d(p ) 


dp>  0, 


-q »dj)  >  0 . 


Then  from  (19)  and  (20), 


A  fo  i  a  d  &  A  „  ^  ~ 

o:L  =  Aa: =  — o  :  o  >  0 

da  2a 


(21) 

(22) 


(23) 


(24) 


rE  dys 


-J  p—!—D  = 
dD 


2(Pl 


2Kr 


>0, 

(25) 


meaning  that  the  kinetic  relations  (19)  and  (20)  are  consistent  with  the  laws  of 
thennodynamics  (23).  Upon  consideration  of  (22),  the  following  constraint  follows  in  terms 
of  the  pressure,  since  (p  >  0  from  (20): 


JE(PlP  -  (l  +  ) (^  /  2  +  K23e  / 3  +  /  4) .  (26) 

Material  parameters  are  selected  so  (26)  is  satisfied  over  the  model's  range  of  applicability. 
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3.  FRAGMENTATION 


Two  approaches  for  predicting  fragment  size  and  velocity  statistics  are  derived  in 
what  follows.  The  first  method,  termed  the  'energetic  approach',  relies  on  a  local  energy 
balance  to  determine  the  mean  size  and  number  of  fragment(s)  associated  with  each  local 
volume  element  (e.g.,  a  finite  element,  computational  cell,  or  particle  in  a  numerical  context). 
The  local  velocity  of  that  volume  element  is  then  assigned  to  all  fragments.  The  second 
method,  termed  the  'statistical  physics  approach’,  is  based  upon  maximization  of  a  global 
statistical  entropy  function  subject  to  constraints  on  conservation  of  mass,  energy,  and 
momentum.  The  premise  of  the  latter  method  is  that  the  most  chaotic  distribution  consistent 
with  the  fundamental  laws  of  mechanics  and  thermodynamics  is  deemed  most  favorable. 
Though  not  pursued  here,  the  latter  stochastic  method  could  be  used  in  conjunction  with  the 
energetic  approach  for  describing  subscale  fragment  distributions  smaller  than  the  grid  size. 

3. 1  Energetic  approach 

Fragmentation  of  the  local  material  element  takes  place  over  a  time  increment 
beginning  when  the  damage  reaches  a  threshold  value,  D  =  DT.  The  energy  released  per  unit 
volume  over  this  time  period  due  to  internal  micro-cracking  is,  from  (15), 


^CLdD  = 
2  Kc 


Kccoc 
2  Kb 


(D-Dt). 


(27) 


During  the  fragmentation  event,  the  material  retains  all  energy  apart  from  eD ,  which  may  be 
stored  or  dissipated  as  heat.  Prior  to  fragmentation,  energy  released  via  r  contributes  to 
dissipation  and  temperature  rise  through  (12);  then,  upon  D  =  DT,  sufficient  gaps  between 

micro-cracks  are  assumed  to  have  developed  such  that  this  energy  is  contributes  to  the 
relative  kinetic  energy  of  fragments.  The  energy  transfer  from  the  bulk  constitutive  model  to 
fragment  expansion  then  follows  as 


MeD=Mplf\D ^  =d'^j(ul+us)/ dt ,  (28) 

where  M  is  the  total  mass  of  the  volume  element,  and  u,  and  us  denote  absolute  energy  of 

the  volume  due  to  contributions  from  the  relative  linear  and  spin  momentum  per  fragment. 
The  progression  of  damage  and  fragmentation  is  illustrated  in  Fig.  3. 


Undamaged 
D  =  0 


Damaged 
0  <  Z>  <  Z)r 


Fragmenting 
Dt<D<  1 


crack  energy  and  friction 
— >  temperature  rise 


crack  energy  release 

— *■  fragment  expansion  KE 


Fig.  3.  Damage  evolution  and  fragmentation  processes. 


For  cube-shaped  fragments  with  linear  edge  dimension  b,  energies  in  (28)  can  be  estimated  as 
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Uj  =  b5  pe2  / 1 6  , 


us  =  b5  pc/)1  / 12 , 


(29) 


where  the  scalar  strain  rate  in  the  fragment  isi  =  V D :  D  ,  with  D  the  symmetric  part  of  L,  and 
(j)  is  the  rate  of  rotation  of  the  fragment  about  its  moment  of  inertia.  Derivation  of  the  first  of 
(29)  follows  that  in  [10].  Combining  (28)  and  (29), 

(^pb3)eD=^(b5pe2/l6  +  b5p</>2n2).  (30) 

Upon  assuming  equal  and  constant  mass  density  among  N  fragments,  replacing  b  with  a 
median  effective  length  b 3  =  ( pN)  1  M  ,  (30)  is  integrated  to  yield 

(31) 

Equation  (31)  produces  the  mean  fragment  dimension  b  .  Mass  conservation  then  provides 
the  local  number  of  fragments  N  =  M  /  (pb3^. 

3.2  Statistical  physics  approach 

Here,  the  entire  fragmenting  body  is  considered  at  once,  i.e.,  globally.  The  mass 
distribution  is  predicted  via  entropy  maximization,  subject  to  the  constraints  that  the  total 
mass  and  total  number  of  fragments  are  known,  following  [11].  The  constraints  are  written 
as 

M  =  pYJbi=pNb\  (32) 

where  the  mass  density  p  is  assumed  the  same  among  all  fragments  and  the  linear  velocity 

does  not  affect  the  mass  distribution.  Quantities  M  and  b  are  determined  from  mass 
conservation  and  a  global  application  of  the  preceding  energetic  analysis,  respectively.  The 
statistical  entropy  associated  with  the  mass  distribution  is  given  by  SM  =  kB  In  P ,  where  kB  is 
Boltzmann’s  constant  and  P  is  the  number  of  conceivable  fragment  arrangements.  The  value 
of  SM  is  maximized  via  maximization  of  the  quantity  [17] 

In  P  =  N  In  N  -  ^  ni  In  nt ,  (33) 

i= 0 

where  nj  is  the  number  of  fragments  of  mass  mi  <  m  <  mi  +  8mi ,  with  8mi  describing  the 

range  of  masses  admitted  in  each  bin  i,  where  j  spans  the  total  number  of  bins.  Constraints 
(32)  are  rewritten  as 

N  =  Yjni ,  M=Yjmi.  (34) 

1=0  i= 0 

Introducing  Lagrange  multipliers  a  and  /3  ,  an  equivalent  function  /(«,  )  is  defined  as 

f  =  \nP  +  a(N-Yjni)  +  j3(M~Yjmi),  (35) 

maximized  by  the  solution  of 
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dni 


(36) 


=  -  (l  +  In  nt )  -  a  -  [5mi  =  0 . 


With  the  notation  change  a  + 1  — »  a  ,  (36)  yields 

ni  =  exp  (-a  -  ft mi ) .  (37) 

Rewriting  (34)  as  integrals  of  continuous  functions  [11]  gives 

oo  oo 

N  =  J ndm  =  J exp(-a  -  f:hn j dm  =  exp(-a)/ /?  ,  (38) 

0  0 

00  00 

M  =  J  mndm  =  J  m  exp  (-a  -  f:hn)dm  =  exp(-a)/ J32 .  (39) 

0  0 

Finally,  the  mass  probability  distribution  function  (37)  becomes 

n  ( mi  )  =  nj  ={N2  / M )  exp  ( -Nmi  / M ) .  (40) 

The  cumulative  probability  distribution  of  fragments  larger  than  m  is  then 

n(m)/N  =  N0  exp(-Ar0m) ,  (41) 


where  N0M  =  N  .  Cumulative  distribution  (41)  is  shown  in  normalized  form  in  Fig.  4(a). 


Fig.  4.  Cumulative  mass  distribution  (a)  and  velocity  probability  distribution  (b),  statistical  physics  model. 


Analogously,  the  measure  of  statistical  entropy  of  the  fragment  velocity  distribution  is 
Sv  =  kB  In  W ,  attaining  its  greatest  value  upon  maximization  of 

In  W  =  N  In  N  -  ^  nj  In  nt ,  (42) 

(=0 

where  nj  is  the  number  of  fragments  with  kinetic  energy  ei<e<ei+8ei,  with  8e/ 

describing  the  range  of  energies  admitted  in  each  bin  i.  Constraints  on  the  fragment  velocity 
distribution  are 

JV  =  £»,>  £  =  1>,,  (43) 

i=0  i= 0 


with  E  the  total  kinetic  energy  of  the  fragment  cloud  that  will  be  determined  later.  Applying 
an  analogous  procedure  as  in  (35)-(37),  the  velocity  distribution  takes  the  form 
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ni  =  exp  (-a- /?e(.) , 


(44) 


where  Lagrange  multipliers  a  and  /?  are  detennined  as  follows.  Let 

j3  =  \/kBT  =  3N/2E,  (45) 

where  T  is  a  thermodynamic  temperature,  via  analogy  with  rigid  molecules  [17].  Assume 
that  in  ballistic  scenarios  the  fragment  velocity  is  coaxial  with  that  of  the  center  of  mass  of 
the  fragment  cloud,  such  that  2 <2.  =  mtvf .  Then  the  velocity  distribution  may  be  expressed  as 

h  ( v, )  =  A  exp  ( -3 Nm^f  /  4A7  ^ ,  (46) 


with  ^4  =  exp(-«)  determined  by  normalization: 

00  00 

1  =  | h{v^dv  =  A J exp(-3 Nmv2  / AE^dvi  =  A^tvE /3mN  , 
0  0 

giving 

n  (v)  =  V3 mN / nE  exp (-3 Nmv2  / 4fs ) . 


(47) 

(48) 


The  form  (48)  differs  from  that  obtained  by  Grady  &  Winfree  [12],  who  assumed  fragments 
may  scatter  randomly  in  all  directions,  yielding  3D  Maxwell-Boltzmann  type  velocity 
statistics.  In  the  ballistic  events  simulated  later,  the  flying  debris  tend  to  follow  a  path  closer 
to  one-dimensional  than  3D  or  radial.  Velocity  probability  distribution  (48)  is  shown  in 
nonnalized  fonn  in  Fig.  4(b). 

The  joint  probability  distribution  of  mass  and  velocity  is  derived  by  combining  (43) 
and  (48): 


p(m,v )  =  h  ( m )  n  ( v) 


3  mN5 
V  xEM2 


exp 


(  N  3  N  2V 

—m 

—  + — V 

(M  4  E  ) 

(49) 


Then, 


J7H  m,  v)  mvdmdv  =  J  mft  ( in )  j  J*  n  (v)  vdv  dm  =  4 EM  /3 


0  0 


(50) 


Identifying  (50)  with  the  linear  momentum  of  the  fragment  distribution,  conservation  of 
linear  momentum  dictates  that  the  energy  partitioned  to  the  velocity  distribution  is 


E  =  3pbim2.  (51) 

Note  that  the  expansion  energy  of  (27)  does  not  contribute  to  the  velocity  distribution  because 
this  energy  is  consumed  by  the  strain  and  rotation  rates  of  the  fragments  as  indicated  in  (30). 

4.  CONCRETE:  PARAMETERS  AND  IMPLEMENTATION 

The  preceding  constitutive  theory  was  applied  to  describe  a  concrete  of  unconfined 
compressive  strength  48  MPa  (7  ksi),  as  studied  previously  in  [1,  4].  The  fully  crushed 
material  is  assumed  to  behave  like  the  aggregate  at  high  pressures,  following  Holmquist  et  al. 
[4].  Constants  Aj ,  K2,  and  A7 ,  describing  the  volumetric  elastic  response  of  the  fully 
crushed  material,  were  determined  by  fitting  the  pressure -volume  response  of  (17)  to  the 
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shock  Hugoniot  data  for  granite  [18],  assuming  volumetric  elastic  deformation  of  the  fonn 
FE=(l-/)l.  Porosity  evolution  parameter  a  was  found  via  linear  interpolation,  i.e., 

a  =  (pLcQ  /(pL- pc).  Thennodynamic  admissibility  of  porosity  evolution,  from  (26),  was 
verified  numerically  for  all  <p>  0  . 

Fracture  toughness  Kc  was  obtained  from  ref.  [19].  Following  Sinha  et  al.  [20], 
Hanchak  et  al.  [1],  and  Hohnquist  et  al.  [4],  failure  occurs  ( D  =  \)  when  the  cumulative 
plastic  strain  A  =  0.0033  under  null  pressure  conditions  p/a0=  0,  giving  £  =  300. 
Similarly,  following  from  data  and  analysis  in  the  above  references,  prescribing  D  -  1  when 
A  =  0.01  under  average  pressure  p  /  a0  =1/6  yields  nD  =  4  . 

Maximum  crack  density  coc  and  threshold  damage  parameter  Dr  were  chosen  based 

on  the  assumption  that  the  typical  fragment  size  is  on  the  order  of  the  minimum  dimension  of 
the  coarse  aggregate  of  the  concrete  microstructure,  here  9.5  mm.  For  cubic  fragments,  this 

implies  an  edge  length  b  of  9.5/  V3  mm.  A  typical  strain  rate  observed  over  the  duration  of 

the  fragmentation  event  was  ,7  =  2(1  0)  /s ,  based  upon  results  of  numerical  analysis  to  be 

discussed  later.  With  these  values  of  b  and  s ,  invoking  (31)  without  spin,  and  assuming 
that  fragmentation  begins  at  the  partially  damaged  state  characterized  by  DT  =  0.5  ,  gives  the 

critical  micro-crack  density  coc  =  1.7(105)  m'1 .  Though  the  model  functions  adequately  for 

the  present  investigation,  more  experiments  and  simulations  are  clearly  needed  to  enable 
unique  selection  of  parameters  describing  energetics  of  fracture  and  fragmentation  over  a 
range  of  loading  rates,  stress  states,  sample  sizes,  and  concrete  microstructures. 

The  deviatoric  stress  rate  expression  was  derived  by  differentiating  (16)  and  assuming 
small  elastic  stretch: 

if  =  2G(1-D)De  +  We6-gWe  -Z)6/(l-D) ,  (52) 

where  DE  is  the  symmetric  and  deviatoric  part  of  LE  and  WE  is  the  elastic  spin  (skew  part 
Le).  The  small  elastic  strain  assumption  is  thought  to  be  justified  for  concrete  deforming  in 
shear  or  tension  as  yielding  and  failure  of  the  material  should  occur  prior  to  the  attaimnent  of 
large  deviatoric  stresses.  At  large  pressures,  the  hydrostatic  response  is  integrated  in  terms  of 
the  modified  total  volumetric  strain  ju ,  following  the  approach  in  [4]: 

P  =  (p~(Pl)/{1  +  <Pl)-  (53) 

At  pressures  beneath  the  locking  pressure,  (17)  is  used,  whereas  a  direct  cubic  fit  to  the 
Hugoniot  data  [18]  is  invoked  at  high  pressures.  Possible  errors  due  to  linearization  of  the 
elastic  strain  propagated  from  (52)  are  thus  avoided  in  the  hydrostatic  response. 

The  energy  balance  (12)  is  exercised,  along  with  thermodynamic  expressions  (22)  and 
(25),  to  update  the  temperature  and  internal  energy  of  the  material.  The  elastic  strain 
quantities  of  (5),  and  hence  the  elastic  deformation  gradient  FE ,  are  needed  for  this  purpose. 
The  latter  is  integrated  numerically  via  [21]: 

F'„,=exp(L'A()F'„  (54) 

where  At  is  the  time  increment  of  integration. 

The  following  methodology  for  addressing  failure  of  the  material  is  applied.  When  an 
integration  point  achieves  a  critical  value  of  damage,  D  -  1 ,  'failure'  occurs.  Failed  material 
supports  no  deviator  stresses  or  tensile  pressure.  Finite  elements  are  converted  into  particle 
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nodes  when  a  scalar  measure  of  strain,  termed  the  ’erosion  strain',  is  attained  [7,  8].  In  the 
calculations  that  follow,  the  erosion  strain  is  chosen  as  0.5,  as  recommended  in  [22]. 

Recall  that  two  methods  have  been  developed  for  predicting  fragment  mass  and 
velocity  distributions.  For  the  approach  of  Section  3.1,  a  typical  fragment  dimension  is 
computed  for  each  converted  particle  using  Eq.  (31),  where  the  strain  rate  s  during  the 
fragmentation  event  is  time-averaged  only  over  integration  cycles  for  which  D  >  Dr.  Each 

fragment  is  assigned  the  velocity  of  its  parent  particle.  For  the  statistical  physics-based 
approach  of  3.2,  the  total  mass  and  mass-weighted  velocity  of  the  particle  cloud  (see  later  Eq. 
(56))  are  extracted  from  the  simulation  output.  The  total  number  of  fragments  follows  from 
assumption  of  a  nominal  fragment  dimension  corresponding  to  the  aggregate  size  as 
discussed  above.  Velocity  and  mass  probability  distributions  are  then  calculated  in  a  post¬ 
processing  step. 

5.  BALLISTIC  SIMULATION 

The  initial  problem  geometry  is  shown  in  Fig.  5,  in  which  a  tungsten  sphere  impacts  a 
concrete  plate  at  a  normal  striking  velocity  of  1120  m/s.  Two  meshes  were  considered,  a 
coarse  grid  with  102144  composite  tetrahedral  elements  [22]  and  a  fine  grid  with  244512 
elements.  The  analysis  was  conducted  over  150/zs.  Contact  between  projectile  and  target 
was  assumed  frictionless  in  the  calculations.  The  alloy  of  the  spherical  projectile  was 
modeled  using  the  plasticity  formulation  of  Johnson  &  Cook  [5]  with  a  Mie-Gruneisen 
equation-of-state  for  hydrostatics  [22].  Failure  of  the  tungsten  material  was  suppressed. 


(a)  (b) 

Fig.  5.  Ballistic  problem:  coarse  mesh  (a)  and  fine  mesh  (b). 


Results  from  the  two  simulations  are  compared  quantitatively  in  Table  1.  The  hole 
diameter  was  estimated  from  the  visible  perforation.  The  crater  diameter  was  estimated  from 
the  region  surrounding  the  perforation  where  the  material  was  fully  damaged,  i.e.,  where 
£>  w  1 .  Presumably  such  fully  damaged  material  would  simply  fall  off  the  target  after  the 
impact  event  due  to  gravimetric  force,  leaving  behind  a  crater.  The  mass  loss  in  the 
simulation  was  found  from 

ML=j  pDdV ,  (55) 

where  the  domain  of  integration  in  (55)  is  the  volume  of  target  material.  Eroded  mass  M 
reported  in  Table  1  was  calculated  as  the  summed  mass  of  all  particle  nodes.  The  average 
velocity  magnitude  of  the  fragment  cloud  was  computed  by 

v=M~lYJmkvk’  (56) 
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with  the  velocity  vk  of  each  particle  k  weighed  by  that  particle's  mass  mk . 


Table  1.  Numerical  results. 


Model, 

coarse 

mesh 

Model, 

line 

mesh 

Penetrator 

VR  [m/s] 

814 

825 

Hole 

diameter  [mm] 

25 

26 

Crater 

diameter  [mm] 

52 

50 

Mass  loss 

M L  [kg] 

0.134 

0.123 

Eroded  mass 

M  [kg] 

0.014 

0.014 

Avg.  fragment 
speed  v[m/s] 

133 

116 

Velocity  [m/s] 
1200  _ 


Fig.  6.  Simulation  at  150  jus.  Particles  scaled  by 
fragment  diameter  and  colored  by  velocity  magnitude. 
Target  elements  colored  by  damage  D. 


The  coarse-meshed  simulation,  150  //s  after  impact,  is  shown  in  Fig.  6. 
Corresponding  results  for  the  fine  mesh  were  similar  and  are  not  shown.  Concrete  particles 
are  scaled  visually  by  the  nominal  fragment  diameter  b  ,  as  computed  from  (32),  assuming 
0  =  0.  Particles  are  colored  by  velocity  magnitude.  The  projectile's  exit  velocity  exceeds 
that  of  most  of  the  particles.  Larger  fragments  are  slower-moving  and  remain  close  to  the 
target,  with  smaller,  faster-moving  fragments  in  some  cases  outpacing  the  projectile. 
Contours  of  damage  D  are  shown  in  the  elements  comprising  target  concrete  material.  The 
most  damaged  concrete  logically  surrounds  the  perforation,  though  the  damage  pattern  is  not 
purely  axi-symmetric. 

Mass  probability  distributions  for  the  fragment  cloud  are  shown  in  Fig.  7(a).  The 
distributions  were  computed  in  two  different  ways:  the  energetic  theory  of  3.1  and  the 
statistical  physics-based  theory  of  3.2.  For  the  former,  the  mass  of  fragments  comprising  a 
particle  k  was  found  from  (3 1)  as 


K2ccoc  (D-  Dt  ) 
8  KEpmk1 


(57) 


Distributions  were  then  constructed  by  grouping  fragments  into  bins  organized  by  mass: 

p(ml  <  m  <  m2 )  =  M~l  ^  mj ,  (58) 

J 


where  j  includes  fragments  with  masses  in  the  range  mx  <  m  <  m, .  Probability  distributions 
for  the  statistical  theory  were  found  from  integrating  Eq.  (49)  as  follows: 
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m- 


(59) 


■y  m2  00 

p(mx  <  m  <  m2 )  =  —  J  J mp{m,\)dvdm 


(  m  1  Z 

f-N  } 

- 1 - 

exp 

- m 

{M  N) 

l  M  J 

Minor  differences  in  mass  distributions  among  simulations  with  coarse  and  fine  meshes  are 
predicted  by  the  energetic  theory,  a  few  percent  at  most  for  any  bin.  Since  an  identical  value 
of  eroded  mass  M  was  found  in  both  simulations  (Table  1),  the  statistical  theory  predictions 
were  identical  for  coarse  and  fine  grids. 

Velocity  probability  distributions  are  shown  in  Fig.  7(b).  For  the  energetic  theory,  the 
magnitude  of  mass-weighted  velocity  probability  for  any  particular  bin  was  found  from 


p(vl<v<v2)  =  M1Y,mk’ 

k 


(60) 


where  k  identifies  a  particle  with  velocity  in  the  range  vx  <  v  <  v2 ,  with  mass  mk .  For  the 

statistical  physics-based  theory,  the  center  of  mass  velocity  v  of  (56)  and  Table  1  was 
substituted  into  (51)  to  compute  the  global  kinetic  energy  E.  The  probability  was  then 
estimated  from  (49)  as 


H 


v2  oo 


V,  <  V  <  V. 


)  =  N  ’JJ  p(m,  v)dmdv  =  V-3/V /  4E  (v/ \Jn / M  +  3Nv' 


74  E 


(61) 


Figure  7(b)  shows  the  predicted  mass  fraction  of  fragments  within  a  particular  velocity  range. 
The  bin  of  lowest  velocities,  with  range  0  <  v  <  100  m/s ,  was  predicted,  for  both  methods 
and  both  meshes,  to  contain  the  largest  mass  fraction  of  fragments. 


■■i  Energetic  theory,  coarse  grid 
Energetic  theory,  fine  grid 
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(a) 

It  is  suggested  that  fragment  kinetic 
energy  may  be  of  importance  with  regard  to 
human  lethality  from  exposure  to  the  flying 
debris.  Kinetic  energy  distributions  are 
shown  in  Fig.  7(c),  where  p(ex<e<e2)  is 
the  mass-weighted  probability  of  a  fragment  i 
having  a  kinetic  energy  et  =  myf  /  2  in  the 
range  ex  <  e  <  e2 .  Distributions  were  found 
with  the  energetic  and  statistical  methods 
analogous  to  those  discussed  above  for  mass 
and  velocity.  Note  that  dependence  on  mesh 
density  is  exacerbated  by  the  squaring  of 
velocity  in  the  kinetic  energy  computation. 


v2  [m/s] 

(b) 


■■  Energetic  theory,  coarse  grid 

^■1  Energetic  theory,  fine  grid 
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_ j _ ft  1  ll  r  " _ l _ 

10  100  1000 


e2[J] 
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Fig.  7.  Predicted  fragment  mass  (a),  velocity  (b), 
and  kinetic  energy  (c)  distributions  at  150  /;s. 
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6.  CONCLUSIONS 

A  constitutive  framework  has  been  formulated  with  the  following  features:  finite 
deformation  kinematics  of  elasticity,  plasticity,  and  pore  collapse;  thermodynamic 
contributions  from  nonlinear  elasticity,  damage,  and  porosity;  and  kinetics  consistent  with 
both  the  energy  balance  and  entropy  inequality.  Two  thermodynamically  motivated  methods 
have  been  derived  to  compute  fragment  size  and  velocity  distributions.  The  theory  was 
implemented  in  a  Lagrangian  finite  element  setting  with  GPA,  and  was  used  to  study  the 
impact  of  a  metal  projectile  upon  a  thin  concrete  barrier.  Residual  projectile  velocities,  hole 
and  crater  dimensions,  and  target  mass  loss  were  comparable  between  two  simulations  of 
different  mesh  density.  Fragment  debris  characteristics  appeared  realistic,  with  the  largest 
particles  usually  flying  slowest.  Fragment  distributions  were  qualitatively  similar  among  all 
model  outputs,  with  the  bin  corresponding  to  the  largest  mass  fraction  of  fragments  having  a 
particular  mass,  velocity,  or  kinetic  energy  range  the  same  regardless  of  method  or  mesh  size. 
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